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^-^ Stability perserving is an important topic in approximation of systems, 

CNJ e.g. model reduction. If the original system is stable, we often want the ap- 

{3JQ proximation to be stable. But even if an algorithm preserves stability the 

^ resulting system could be unstable in practice because of round-off errors. 

■^ Our approach is approximating this unstable reduced system by a stable sys- 

T— I tem. More precisely, we consider the following problem. Given an unstable 

^^ linear time-invariant continuous-time descriptor system with transfer func- 

rj tion G, find a stable one whose transfer function is the best approximation 

Q of G in the spaces RH2 and RHoo, respectively. Explicit optimal solutions 

^ are presented under consideration of numerical issues. 

-(— > 

B 

> 



> 



1 Introduction 



We deal with linear time-invariant descriptor systems described by differential algebraic 
(N| equations (see e.g. [KM06], [Dai89] for more details) 

in 

0\ E±(t) = Ax(t) + Bu(t) 

cn yit) = Cx{t) + Duit) 

o 

r^ with matrices {E,A,B,C,D) € M"''" x M"^" x M"^™ x RP^" x W"^ and the corre- 

sponding transfer function s i— )• C {sE — A)^^ B + D . In many applications such as model 
reduction (see [ASGOl] for an overview an further references) or system identification 

^ (see [Lju99]) one wishes to approximate systems such that the transfer functions of the 

3 original and the approximated system are as close as possible. An important property 

of an approximation technique is preservation of stability (cf. e.g. balanced truncation 
[ASGOl], hankel norm approximation [Glo84]). However, not every method has this 
property in general (Arnoldi method, Lanczos method [ASGOl]). Moreover, a stable 
approximation can be unstable in practice [BFF98]. This happens in computer imple- 
mentations because of round-off errors. One way out is restarting the algorithm with 



changed parameters, e.g. interpolation points in Lanczos method. The disadvantage is at 
least doubled computational complexity. Another approach is to modify the algorithm 
directly (e.g. partial Pade-via-Lanczos method [BFOl]). 

In this paper our approach for stability perserving is to consider the computed unstable 
approximation of the original system. Given is an unstable descriptor system S (1.1) 
with transfer function G. Find a stable descriptor system whose transfer function is the 
best approximation of G in the spaces RL2 and -R-Lqo of real rational functions on the 
imaginary axis iM in the Lebesgue spaces L2 and Lqo, respectively. We call this system 
an optimal RH2-a,ppiox.iraation and optimal i?ffoo-approximation of S, respectively. 

For special cases this problem was already solved. The characterization of the unique 
optimal i?/:f2-approximation of a standard system (i.e. £^ = / in (1.1)) with D = is 
a simple consequence of the Paley- Wiener theorem. The set of all suboptimal RH^o- 
approximations of standard systems was determined in [GGLD90] with j-spectral fac- 
torizations to solve model matching problems. However, no optimal ones were explicitly 
given. An explicit representation of a nonunique optimal i?i?oo -approximation of mini- 
mal antistable standard systems was presented in the proof of [Glo84, Theorem 6.1]. This 
theorem was only used as an auxiliary result for a more general problem, the optimal 
hankel norm approximation. That representation requires the possibly ill-conditioned 
balanced minimal realization of G. This numerical issue was discussed in [SCL90]. 

Our approach for approximating an unstable descriptor system S in a numerically reli- 
able way is the following. First, we decompose S into a stable and antistable part by 
combining [GL96, Theorem 7.7.2] and [KD92, Section 4.1]. We show that it is sufficient 
to replace the antistable part by its approximation to get an approximation of S. This 
was e.g. proposed in [Fra87, Section 8.3] for standard systems and w.r.t. RL^o- The 
stable part is the unique optimal i?i:f2-approximation of S. Our representation of the 
optimal iii^oo-approximation only requires at most one singular value decomposition. 
For discret-time standard systems [MarOO] presented an approximation method similar 
to ours. However, in contrast to our approach they used a balanced minimal realization 
of the antistable part. 

The paper is organized as follows. In Section 2 we summarize some results regarding the 
theory of matrices and matrix pairs. We introduce the setting and recall some system 
transformations, e.g. balanced realization. Section 3 states and solves both approxima- 
tion problems. Section 4 presents an algorithm solving our problem and some numerical 
examples. 

2 Preliminaries 

We define C>o := {s G C; Re(s) > 0} and analogously C>o, C<o, C<o. The set iR is the 
imaginary axis. For a function / we denote its domain by T)[f) and for M C P(/) the 
image of M under / by /[M]. The kernel of a matrix A G C"-x™ ig denoted by ker(^), the 
adjoint by A* , the spectrum by cf{A) and the Frobenius norm by ||^||f = Y^trace(A*^). 



The matrix /„ G M"^" is the identity matrix and Onxm G M"^"^ the zero matrix. The 
vector Cj G M" is the zth unit vector. 

Definition 2.1. We call {E,A) £ M"^" x M"^" a matrix pair. A scalar A G C is called 
a (generalized) eigenvalue of (ii^, A), if det(A-E' — ^) = 0. If there exists v ^ such that 
Ev = and Av ^ 0, then oo is called a (generalized) eigenvalue of {E, A). We denote the 
set of all eigenvalues of {E, A) by cr{E, A) C CU{oo} and the resolvent set by p{E, A) := 
C\a{E,A). If p{E,A) = 0, then {E,A) is called singular otherwise regular. We denote 
the set of all regular matrix pairs by P„ := {{E,A) G M"^" x M"^"; {E,A) regular}. 

Remark 2.2. If the matrix pair is regular, then the set of eigenvalues is finite. Otherwise 
every A G C is an eigenvalue. Details about the generalized eigenvalue problem can be 
found in [Ste72]. 

To prepare decompositions of systems, we combine the results in [GL96, Theorem 7.7.2] 
and [KD92, Section 4.1] to get the following theorem. 

Theorem 2.3. Let {E,A) G P^ and Mi,M2 C C U {oo} be two disjoint sets such that 
a{E,A) = Ml U M2. Then there exist orthogonal U,V e M"""" such that 

a(^i,yli) = Ml and cr(^3, A3) = M2. Moreover, there exist R,L e W^^p such that 

AiR - LA3 = -A2 

(2.2) 
E\R — LiE^ — — E2 ■ 

The matrices P := {^^-j^) U, Q:=V {^f) satisfy 

'El 0\ fAi 



iPEQ,PAQ)-,,^ ^^^,^^ ^^ 

Remark 2.4. We get (2.1) e.g. by applying the generalized Schur decomposition to {E, A). 
Equation (2.2) is called generalized (or coupled) Sylvester equation (see [KD92]). 

We recall the following known result which can be used to test regularity of matrices. 

Theorem 2.5 ([ZDG96, Section 2.3]). Let A G C"^", B G C"^"', C G C"'^", D G 
(j~.mxm and M := {^ f)) . If A is regular, then we have det{M) = det{D-CA-^B)det{A). 
If D is regular, then we have det(M) = det(A — BD~^C)det{D). The matrices D — 
CA^^B and A—BD~^C are called the Schur complement of A and D in M , respectively. 

Definition 2.6. For n,p,m G N we define the sets of systems (compare (1-1)) 
S^p^^ := p„ X M"^" X RP^" X RP^™, 
Sn,p,m := {{E, A, B, C, D) G S„,p,^; iR C p{E, A)} , 
Sn,p,m := {{E, A, B, C, D) G S„,p,^; C>o C p{E, A)} , 
Sn,p,m := {{E,A,B,C,D) G Sn,p,m\ C<o C p{E,A), E regular} , 



where S^^p^m ^^'^ ^n,p,m ^^^ the set of the stable and antistable systems, respectively. 
We call a system S = {E,A,B,C,D) £ ^n,p,m a standard system if E = I and a 
descriptor system otherwise. For Si = {Ei, Ai, Bi,Ci, Di) G Sni,p,m, i £ {1)2}, and 
regular P,Q€z M"i^"i we define the operations 

P • 5i • Q := (P^iQ, PAiQ, PBi,CiQ, Di) , 

P ■ Si := P ■ Si ■ Ini , Si ■ Q := In^ ■ Si ■ Q . 

Some of the following definitions and theorems for standard systems can be given for 
descriptor systems, too. However, we restricted ourselves to the specializations we need 
here. 

Definition 2.7 ([Dai89, Definition 2-6.1]). Let S = {E,A,B,C,D) € Sn,p,m- Then 

g{S) : p{E, A)3s^^ C{sE - A)-^B + De W"" 

is called the transfer function of S. We denote the set of all p x m transfer functions by 
Tp,m := {GiS); S G UneN^fi,p,m}- A system S G Sn,p,m is called a realization of G{S), if 
g{S){s) = g{S){s) for aU s G V{g{S)) n V{g{S)). We denote the set of all realizations 
of g{S) by E{g{S)). A realization Si G Sni,p,m of g{S) is said to be minimal, if every 
other realization 52 G Sn2,p,m of ^(5) fulfils n2 > ni. 

Remark 2.8. Transfer functions are meromorphic. Obviously, ^(5i©S2) = ^(5'i)+^(5'2) 
andg(P-5-Q) = g{S) hold. 

We want to transform a system 5" = {E,A,B,C,D) into S* = {E,A,B,C,D) with 
^(-S*) = ^(5"), in particular a(E,A) = a{E,A), such that E,A,B,C,D have a special 
structure, e.g. to decompose systems into a stable and antistable part. Note that all 
matrices stay real under the presented transformations. 

Definition 2.9 ([Dai89, Definition 1-3.1]). Two systems Si = {Ei,Ai,Bi,Ci,Di) G 
En,p,m, i G {1)2}, are called restricted system equivalent (for short r.s.e.. Si ~ 5*2), if 
there exist regular P, Q G M"""" such that P ■ Si ■ Q = 82. 

Lemma 2.10. For every S = {E, A, B, C, D) G Sn,p,m there exist 5^ ~ 5, fc G {0, . . . , n} 
and i^ G N such that 

and N" = 0. Moreover, V{g{S)) = p{J) = p{E,A) holds and 



u-l 



Vs G V{g{S)) : g{S){s) = Cj{sl - J)-^Bj + D-Y, s'CnN'Bn ■ (2.3) 



i=0 



Proof. We apply Theorem 2.3 to iE,A) with Mi = a{E,A)\{oo}, M2 = {00} and 



get regular Ei,A-i,P,Q as in the theorem. Choose P := ( ^1 M P, Q : 



I 



I 
A-^ 



Q 



and set S^u := P ■ S ■ Q. The matrix A^ := £"3 is nilpotent because o'{N, I) = {00} 
and thus c^N) = {0}. Equation (2.3) follows from G{S) = G{Sw) (Remark 2.8) and 



Cn{sN - I)-^B 



N 



X^iLo s^CnN''-Bn (Neumann series). 



n 



Theorem 2.11 ([ZDG96, Theorem 3.10, Theorem 3.17]). Let S = {I,A,B,C,D) £ 



^n,p,m- -I hen 



S 



( (Au Ai3 \ 

A21 A22 A23 A24 

A33 

V V A43 AuJ 



B2 




A 



(Ci C3 0),D 



J 



where An e ]R"'^"», rij e {0, ... ,n} fori E {1,2,3,4} and Sm ■= {I,Aii,Bi,Ci,D) is a 
minimal realization ofQ{S). This decomposition is known as the Kalman decomposition. 
All minimal realizations ofQ{S) are r.s.e. to Sm- 

Now we introduce the controllability and observability Gramians. Later these matrices 
will play a central role in our problem of computing an optimal Pii^oo-approximation. 
The next theorem follows easily from [Pen98, Corollary 2]. 



Theorem 2.12 ([Pen98], Lyapunov equation). Let S = {E,A,B,C,D) G S'^^^. 
there exist Xc, Xq G M"^" such that 

AXcE^ + EXcA'^ + BB'^ = , 
A^XoE + E'^XoA + C^C7 = . 



Then 



(2.4) 



They are unique and symmetric. We denote them by Cs '.= Xc ('controllability Gramian 
of system S) and Os '.= Xq ('observability Gramian j. 



Lemma 2.13. Let S = {E, A, B, C, D) G ^n,p,m <^'^d E he regular. Let P,Q e] 
regular. Then G^s = Q^p.sqQ^ and Ds = P~^ OpsqP hold. 

Proof. We prove this only for <ts since Ds is analogously. By definition of Cp.g., 

= PAQ^Lps-qQ^E'^P^ + PEQ(Lp.s.qQ^ A^ P^ + PBB^P^ , 

= p(A{Q<tp.s.QQ'^)E^ + E{Q^p.s.qQ'^)A^ + BbAp^ 



be 



hold. Due to the regularity of P and Q and the uniqueness of the solution of (2.4) 
(Theorem 2.12) the assertion is proved. D 



Theorem 2.14 ([Glo84, Theorem 4.3]). Let S = {I,A,B,C,D) £ S^.p_„. Then there 
exist Sh ^ S, r G N, h £ M>o such that 



'" ^ -hirj ' ''" V -hir 



If Si = S2 are diagonal, then Sb is known as the balanced realization of G{S). 

To check whether A € C is an eigenvalue of {E,A), we need the following easy general- 
ization of a known result for standard systems (e.g. [Mac04, Theorem 5.6.1 c]). 

Lemma 2.15. Let G € Tp^rn have a realization S = {E, A, B, C, D) € ^n,p,'m- Then for 
every finite A G cy{E, A) at least one of the following three statements holds: 

(z) X is a pole of G, (ii) rank (^XE — A B) < n, {Hi) rank I ] < n . 

Proof. Let S^^, ~ 5 be as in Lemma 2.10. We set S,j := {I, J, Bj,Cj,D){£ Sfc,p,m)- By 
Lemma 2.10 we know that A G a{E,A)\{oo} <^ A G cr{J). By [Mac04, Theorem 5.6.1 
c] one of the following three statements holds: 

{i') A is a pole of G{Sj), {ii') rank (A/ — J Sj) < k, {Hi') rank ( j < A; . 

Equation (2.3) shows that {i) <^ {i'). By the regularity of AA^ — I we have 

rank (A/ - J Bj) + {n - k) = rank ( •^ j = rank {XE - A B) 



and analogously rank ( ^^ \ + {n — k) = rank {^^^c )■ Thus {ii) <^ {ii') and {Hi) <^ 
{Hi'). This finishes the proof. D 

Definition 2.16 ([ZDG96, Section 4.3]). We approximate unstable systems with respect 
to the following spaces and corresponding norms. Let p, tti G N. 

. RL^^"^ := {G G Tp,„; iM C V{G), \\G\\^ < 00}, ||G||^ := sup ||G(i^)||2; 



wei 



RH^™' := < G G RL^"^; G analytic on C>o h 



. i?Lf ™ := {G G Tp,„; iM C P(G), ||G||2 < 00}, ||G||2 := J^ f \\G{ico)\\ldij; 

y -00 

• RH2^™' := < G G RL^^; G analytic on C>o >■ 
We omit the indices, if there is no risk of confusion. 



Remark 2.17. These spaces are often introduced as a subset of socalled real rational 
functions, i.e. entrywise a quotient of polynomials with real coefficients. Actually, every 
real rational function is a transfer function (see [Dai89, Theorem 2-6.3]) and vice versa, 
since G{s)ij = eJC{sE — A)~^Bej + Dij is the Schur complement of sE — A in M{s) := 

( e'^c D— I ^'^'^ ^^^^ G{s)ij = det{G{s)ij) = j(f^.A^_j/\ by Theorem 2.5. Let F and G 
be elements of one of these spaces with S{F) = S(G), e.g. s i— ;• -^^ and s i— ;■ . _^w^_2n ■ 
By definition they are equal on the cofinite set 'D(-P) nX>(G). Therefore we identify them 
as the same function. With this and by the identity theorem for analytic functions the 
spaces are indeed normed spaces. 



Theorem 2.18 ([Fra87, Section 2.3]). Let G_ G RL2 be analytic on C<o and G+ G RH2. 

|2 _ wn I|2 I wn I|2 

I2 — ll<-^-|l2 + ll<-^+ll2 



Then ||G_ +G+\\l = \\G-\\l + \\G+\\l holds. 



The following known result characterizes these spaces in terms of realizations of transfer 
functions. 

Theorem 2.19 ([Fra87, Section 2.3]). Let G G Tp,^ with iR C V{G). Then the following 
hold: 

(a) GgRLoo ^3iI,A,B,G,D)£SiG), 

(b) G G RH^ ^ 3 (/, A, B, G, D) G S(G) n S+p_^ , 
(c)G€RL2 ^3{I,A,B,G,0)€S{G), 

(d) G G RH2 ^ 3 (/, A, B, G, 0) G S(G) n S+p^^ . 

3 Optimal stable approximations 

3.1 Problem statement 

In this paper, we deal with the following problem. 



Approximation-Problem (APg). Let q G {2, 00}. For S G S^„m fi''^d ^ ^ 
I Js^M S^ such that 

'^n€Pi n,p,m 

\\G{S)-g{S)\U = ^ inf^ \\g{S)-g{S)\u. 



In the next subsections we show that this problem is solvable. We also present an explicit 
solution in Theorem 3.2 for q = 2 and Theorem 3.16 for q = 00. 

We do not consider systems S = {E, A, B, C, D) G Sn,p,m where {E, A) has imaginary 
eigenvalues. If there is no imaginary pole of g{S), then there exists a realization 5*2 G 
^n,p,m ofg{S) and we solve (APg) for 82- Otherwise, for every S G UfieN^^nm ^^ have 
life's) — ^(S')llq = cxD and hence (APg) is not solvable. 



First, we show that solving (APg) for a descriptor system we can equivalently solve 
(APq) for its antistable part. This was e.g. proposed in [Fra87, Section 8.3] but only for 
standard systems and g = oo. We need Theorem 3.1 to use some results in [Glo84] and 
[SCL90] which are only applicable for antistable standard systems. 

Theorem 3.1. Let S € S^_p^m- Then there exist 5+ = (£'+, A-|-, i?+, C+, D) G '^nj^^p^m 
and S^ = {E^, A^, B^,C-,0) G ^n-,p,m such that 5 ~ 5+ © S'_. Let q £ {2,oo} and 
7 > 0. Then the following two are equivalent: 

(i) 35 G U S±^_^: \\g{S)-g{S% < 7, (ii) 35 G U Sj^^^: ||g(5_) - g(5)||, < 7 • 

nSN neN 

If S satisfies {ii), then S (B S^ satisfies (i). 

Proof. Theorem 2.3 assures the existence of such systems Sa. and 5_. Let 5 G S~ 

■' J -r n,p,m 

satisfy (i). Then we have 

7 > 11^(5) - ^(5)11, = ||e(5_) - {g{s) - g{s+m, . 

Thus S := S (B (£'+,^+,-6+, — C+, — D) satisfies {ii) (compare Remark 2.8). The fact 
that 5 G UnsN^j^^pm follows easily from the definition of ©. Let 5 G S~„^ fulfils (ii). 
Then analogously S Q) S-\- satisfies (i). D 

3.2 Optimal i?i?2-appi'Oximations 

Theorem 3.2. Let S G ^npm- ^^ "^^e the notations of Theorem 3.1. Then 5+ solves 
(AP2). More precisely, we have 

inf , \\g{s) - g{s)h = ms) - g{s+)h = ms-)h ■ 



S^Unem^. 



n,p,m 



The solution is unique in the following sense: If S2 is another solution, then it is also a 
realization ofg{S+). 

Proof Let 5 G UteN^Jp,™ ^ith g{S) - g{S) G RL2. Since g{S-) G RL2 (Theo- 
rem 2.19) we have ^(5+) — g{S) G RH2- By Theorem 2.18 we obtain 

\\gis)-g{s)f2 = ms-)+g{s+)-g{s)g = ms^)f2+ms+)-g{s)g > ||a(5_)||i . 

The lower bound is attained, if and only if 5 G S(^(5+)), e.g. if 5 = 5+. D 

Remark 3.3. In the proof of Theorem 3.2 it can be seen that the function ^(5+) is 
even the best approximation in the Hardy space H2. The reason is that ^(5_) is also 
orthogonal to H2 (see the more general statement of Theorem 2.18 in [Fra87, Section 
2.3]). 



3.3 Optimal i^i7oo-app>'oximations with balanced realization 

For a minimal antistable standard system 5 the problem (APqo) was already solved in 
[Glo84]. Combining [Glo84, Theorem 6.1] and its proof based on [Glo84, Theorem 6.3], 
we obtain the following result. We reformulate the statements in terms of realizations 
instead of transfer functions. Note also the interchange of the stable and antistable 
system which is easy to show. 

Theorem 3.4 ([Glo84]). Let S G S^^p^m ^^ ^ minimal realization such that 



-.,V' ^^=U -.,/J' -i^-(^^^^)' 

where a\ := y^max a {^s^s) (compare Theorem 2.14-)- Then there exists U G RP^*" with 
B2 = -CjU. Set 

r :=SiS2-a?I, 

I := T-^{afAj^ + Sa^nSi + aiCjUBj) , C := CiSi - aiUBj, 

B ■.= r-^{^2Bi-aiCjU), D:=D + aiU. 

Then S := {I, A, B, C, D) G '^n-r,p,m holds and 

^ ini^ \\Q{s)-g(s)u = \\g{s)-g(s)u = cji. 

Remark 3.5. The value o"! is the greatest hankel singular value of Q{S) (see [Fra87, 
Section 5.1] for further details). As stated in [Glo84, Theorem 6.1], the function Q{S) is 
even the best approximation in the more general Hardy space H^. 

Because of Theorem 3.1 we are now able to solve (APqo), if 5* is a descriptor system. The 
disadvantage of algorithms based on Theorem 3.4 is that we need a balanced minimal 
realization of Q{S^) to compute an optimal solution. 

3.4 Optimal i^i^oo-approximations without balanced realization 

Some results of [Glo84], concerning optimal hankel norm approximation, were generalized 
in [SCL90] to nonminimal stable standard systems. We use [SCL90, Theorem 1] and 
its proof to solve (APqo) for general 5 G '^n,p.m- Among some modifications, one of our 
main tasks is to prove that the optimal solution lies indeed in S^^p^^i since [SCL90] only 
considers the poles of the resulting transfer function. 

First, we need to compute the greatest hankel singular value ai of Q{S) with the formula 
given in Theorem 3.4 which is restricted to minimal standard systems. However, we 



avoid the determination of a minimal realization. To the author's knowledge, it was 
never shown that the nonzero hankel singular values of Q{S) equal the roots of the 
nonzero eigenvalues of ^s^s independently of the realization S of G{S). This equality 
was only proved for minimal realizations (cf. e.g. [Fra87, Section 5.1, Theorem 3]). 

Lemma 3.6. Let S = {I, A, B, C, D) G S^ ^ be a minimal realization of G G Tp^m o,nd 
let S = {E, A, B, C, D) E Sr be a realization of G. Then 

In particular, the equation af = maxcr(C^g£' Dgi?) holds for o-^^ S G S-^^H S{G). 

"■^^ such 



n,p,'m 



Proof. We set Sj := E ^ ■ S. By Theorem 2.11 there exists a regular T G 
that 



/ /^ 



11 



A 



13 



o\ 



St ■.= T-Si- T-^ 



L 



A21 A22 A23 A24 
A33 
V V A43 ^44/ 



/51\ 

B2 





\ 



,(Ci C3 0),D 



I 



where S\ := {I,Aii,Bi,Ci,D) is a minimal realization of G. The condition a(E,A) 
a{E-^A) C C>o implies 



a 



11 



A 

A21 A 



22 



c 



'>0 



a 



An A 

A 



13 



33 



C 



'>0 



Xi X2 



By Theorem 2.12 there exist unique solutions I ^t ^^ ) ^^"^ 

T 



ct(Aii)^C>o 



All 

^21 ^22 

All ^13 

A33 



) C§ S) H.1 ^.:) (t I) Hi)w-.^)-> 



11 1^2 



+ 



Yi Y2 

y/ y4 



All ^13 

A33 



+ 



h]{Ci C3) 



0. 



Thus X 



GTs, , yi = Ds, , 



(Xi X2 0\ 

X^ X4 



V oy 



and Dc 



/i^i y2 o\ 
0000 

y2^ y4 

V oy 



By Theorem 2.11 and Lemma 2.13 there exists a regular P € M"^" with £5 



p-l^ p~T 



10 



Ds^= P^ OsP- We conclude 



a{iS:sE'^DsE)\{0} = a{€s,Os,)\{0} = aiT^Cs^T'^T-^ Os^T-')\{0} 

fXiYi X1Y2 0\ 



a{(ts^Ds^)\{0} = a 



\{0} 



XjYi X^Y2 

V 0/ 

a{XiYi)\{0} = a(e:s,DsJ\{o} = a(e:sDs)\{o} . 



D 



The next theorem is one of our main results. Together with Theorem 3.1 it gives sub- 
optimal solutions of (APoo) and, with an additional regularity-condition, an optimal 
solution. Theorem 3.14 solves (APqo), if the regularity-condition is not fulfilled. 



Theorem 3.7. Let S = {E, A, B, C, D) G S„ ^. Choose 7 > ai and set 



n,.., := D,E(t,E^ - -f^I 



-S,7 



£s,j '■— E 7^5,7 ! 
>lo.^ :^ —A l\.Q., — C, 



^S,7 •- 



-s,-, 



S,7 ! 



Cs,7 



E ' DsB, 
C^,E^. 



If{£s,-y,-^s,-y) is regular, e.g. if'j > ai, thenSg,^ ■= {£s,-y,-As,^,Bs,^,Cs,-y,D) e ^t,p,m o-^d 
Ss~, fulfils 

cTi<\\G{s)-g{SsM^<i. 

Proof. First, let E = I. In this and the following proofs we use the equation 

As,, = -A" lis,, - C^Cs,, = -TZs.^A'' - Bs.^B" = -{" A" + QsA^s. 
which holds because of 

-A^Tls,, - C^Cs,, = I'^A^ - A'Os^s - C^C^s = fA^ + QsA^ 

= -i^A^ - Qs'^sA^ - OsBB'^ = -TZs_^A'^ - B,-,B'^. 



(3.1) 



-S,7^ 



^5,7^ 



Thus As,-, is indeed the same as in [SCL90, Theorem 1]. To apply [SCL90, Theorem 1], we 
first construct an antistable system 5*2 with a hankel singular value greater than 7. This 
method was proposed in [Glo84, Remark 8.4]. Let 71 > 7. We set Si := (1, 571, 71, 71, 0) 
and 

"^'^=^'(0 hJ'VO tJ'VO tJ'VO o))^^n+l,p+l,n.+l- 

Then we have 

-71/ V ~7i. 

Thus 71 is the greatest hankel singular value of G{S2) and cri < 7 < 71. If we apply 
[SCL90, Theorem 1] to 5*2 with iC = G RHoo, then the transfer function of 



'•Sa 







'S2,7 



£. 



S.l 







Yi - J 



2 ' 



A 







S,7 ^ \ { '-'3,1 

^(7^71 +7?)/ I 



Bs 





-7? 



C 



3,1 







-7? 



D 
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is a suboptimal solution, i.e. ||t/(S'2) — ^(532,7) || < 7- Note that [SCL90, Theorem 1] 
does not state, that for every K G RHoo with ||i^||oo ^ 1 the resulting transfer function 
is a suboptimal solution. However, this was shown in the proof (see also [Glo84, Remark 
8.4]). With 

^(^2) - V G{S,)) ' ^^^'--^ - \ G{Ss,, 

we conclude 

\\g{s) - g{Ss,,)\\oo < 11^(^2) - Q{Ss,Moo ^ 7 • 

Now we show that Ss,-i G Sn,p,m- [SCL90, Theorem 1] states that QiSs^n) ^^s at most 

one pole A in C>o counting multiplicities. The pole of t/(5s^,^) is A = 2^^ ^ 3 > with 

7^ — 7^ > 0, 7^71 + 73* > 0. Thus A = A and ^(55,7) has no poles in C>o. 

Let A G C>o. We use Lemma 2.15 to show that A ^ a{£s,'^,^s.'y)- Let u G M" satisfy 
{X£s,-y — •^s,-^)v = and Cs,-,v = 0. Hence with As..y = —A^Ss,'^/ — C~^Cs,^ we get 
(A/ — {—A^))£s^-yV = 0. Because of a{—A^) C C<o we conclude £3,-1'^ = 0. Thus 

V s G C : {s£s,-, - As,-,)v = {s£s,-, + A^ £3^-, + C'^ Cs,^)v = . 

Since the matrix pair (£^5_^,^g_^) is regular we have t; = 0. 

Now let V G M" satisfy {\£j^ - Al^)v = and BJ^v = 0. Hence {XI + A)£j^v = 
^l, = -A£l,-BBl^: 



follows with Al = —A£'J —BBj. Again we deduce v = 0. Therefore A ^ a{£s,j,As,j) 



and thus ^s,^, G 'S>n,p,m follows. 

Finally we consider general regular E G M"^". We already proved that the assertion 
holds for 5^-1 g With (^^-1.5 = C^s and D^-i.g = E~^DsE (Lemma 2.13) we conclude 



£^-i.s„=E'^OsE<ts-7^I, 


B^-i.s,,=E^OsEE-'B, 


^.-.s,7 = -A'E-^£,-.,,,^ - 


-CCg-i.<j_^, Cj,-i.s^ =C(ts ■ 


! see that Ss,^ = S^-i,g^ • E~^ . 


Thus the assertion also holds for Ss,-y ~ '^^-i.^^ 



If 7 > o"i, the matrix pair {£s,^,As,j) is regular, since TZs,-y, E and thus £s,^ are regular. 

D 

We formulate some equivalent conditions for the regularity-condition of [£s,„j^,As,aj^) in 
Theorem 3.7. 

Corollary 3.8. The following three statements are equivalent. 

(i) {£s,^-^,As,a:,) is regular, (ii) a {£s,a:,,As.ai) '^ C<:oU {00}, (Hi) As. ^^ is regular . 

Proof. The implication "(i) =^ (ii)" follows from Theorem 3.7. The statement (ii) implies 
^ a[£s,„j^,As,ai) and thus (0 • £3.^1 + -^s.^i) is regular. Finally, (iii) implies ^ 
a[£s,ai,As,^J, therefore p[£s,aj^,As,^^) / and thus (i) holds. D 
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To get optimal solutions if [£s.a^,'As,ai) is singular, we analyze Ss,j in Theorem 3.7 for 
7 — ^ ai+. For the proofs we use the realization 5^ ~ S* of G{S) as in Theorem 2.14 
with h = (Ti because we will take advantage of the special structure of Ss^,a^- First we 
investigate the relation between Ss,-^ and Ss^^,-^- 



Lemma 3.9. Let S = {E,A,B,C,D) G S'^p^^. Let P,Q G 
Sp.s.Q,-y = Q^ • Ss,., ■ P'^ holds for all 7 G M. 



he regular. Then 



Proof. We set S := P-S-Q. By Lemma 2.13 the equations (Tg = Q(tsQ'^ , Dg = P^OgP 
hold. This implies 

Bs,-, = E^DsB = E^P^DsPB = Q-^Bs,, , 

Cs,, = C(tsE^ = CQ^sQ^E^ = C§,,P-T, 

£s,-y = E'^OsE^sE^ - 7^^^ = E'^P^DsPEQttsQ^E^ - 7^^^ ^ Q-^g^^p-^^ 

As,, = -A^E-~^£s,, + C^Cs,, = -A^E-^Q-'^£s,,P~'^ + C^Cs,,P-^ 

= Q-^As„P~^. D 



Lemma 3.10. Let Sb = {L,Ab,Bb,Cb,D) ^ S be as in Theorem 2.14 with h = o"i. Let 
Ti G M"^" he regular such that 



Ti-S- E-'T^' = Si 



I 0\ (All A 

/J'U21 



tO-dO-e^' ^^).°) 



Then the following three statements hold: 

(a) B2B2 = C2 C2, 

(b) ker(52T) = ker(C2), 

(c) (fs,.i, A,.J singular ^ ker(C72) / {0} <^ }^ei{Bj) / {0}. 

Proof. We get B2B2 = C2C2 by adding the equations —(J1A22 — (y\A22 + B2B2 = 
and (T1A22 + o'i^22 — CJC2 = resulting from (2.4). Therefore 



B^v 



4^ = v'B2BJv = v^CjC2V 



V G ker(B2 ) <^ ^2 u = 4^ 

<^ IIC2UII2 = <^ u G ker(C2). 

By Lemma 3.9 and with F := S2S1 — a\L we have 

Ss„., =T,-^E-^-Ss,.,-T^ 

_((T 0\ (-Aj.V-CjCi^i aiCjC2\ ( ^2Bi\ , . 

"VVO Q)\-Aj2T-C:iCi^i aiC:^C2)\-aiB2)'^^'^' -^1^2), 

Therefore {£s,a^,As,a^) is singular, if and only if (^£s^^,a-^,As^,cT^) is singular. 
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Let ker(C2) 7^ {0}. Then there exists v ^ with C2V = 0, i.e. for all s G C 

/0\ _ (sT + Aj,T + C^CiSi -aiC7C2\ /0\ _ /0\ 

[sts,,., - As,,.,) ^^j - ^ ^T r + cjci^, -aiCjC2) [v) - [0) ■ 

That means [£si^,a-,,Asf^,aj) is singular. 

Let ker(C2) = {0}. Then CJC2 is regular. Let s G C. The matrix sSs^^^i — -^Sj.ai is 
singular, if and only if the Schur complement 

(si + (aJ, + cJCi^iT-' - cJC2icJC2)-\Aj2T + cJCi^i)r-Ay =■ {si - M)r 

of sSs^^rj-, —As,,ai in —<TiCjC2 is singular (Theorem 2.5). This is equivalent to s G a{M). 
Hence (fs^^j,^^^^^) has a finite number of eigenvalues and consequently is regular. D 

Similar to the proof of [SCL90, Theorem 1] we want to eliminate the singular part of 
the system Sg^^i and show that the resulting system is an optimal solution. 

Lemma 3.11. Let (iSs_^^,^s,<ti) be singular. Then there exists a regular T G M"^" such 
that _ _ _ 

-T c rp-i _ ( ( Ell 0\ (All 0\ ( Bi 



where {En, An) is regular. 

Proof. Our proof is based on the proof of [SCL90, Theorem 1]. Let Sb, Ti and T be 
as in Lemma 3.10. By Lemma 3.10 we have d :=dim(ker(C2)) > 0. Therefore there 
exists an orthogonal W G W^'' such that {62 Opxd) = C2W~^ and ker(C'2) = {0} 
(e.g. singular value decomposition). We define {BJ BJ) := BJW~^ with appropriate 
dimensions. By Lemma 3.10 we have BJ = 0. Now we set T2 := (qw')) conclude 
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Bi 

'Jdxm 



A{2r — 6*2 CiSi = aiB2B-l from (3.1) and compute 

\ Odxm / 

Cs„^,T^^ = CbDs.T^^ = (CiSi -aiC2W^) 

= (CiSi —aiC2 Opxd) ='■ [Ci Opxd 



An 
Odxd 



T2^S„.,T2 -T2\^ ^^^^^T aiC2^C2J^2 

-AT^r - C^CiSi aiC^Cs 

£71^2^7 aiCjC2 

Odxdy 



Since C2C2 is regular, there exists s G C such that the Schur complement 
'si + (aJ^ + CTc7iSir-i - aiCjC2{C^C2)-^B2BjT-^\) T 



of S-Eii — All in —(T1C2C2 is regular. Theorem 2.5 states that S-En — An is regular. 
Thus (£^11,^11) is regular. With T := T2T^ the proof is finished. D 

Lemma 3.12. We set Si ■.= {Eii,Aii,Bi,Ci,D). For every uj £R with iu £ p{Eii, An) 
the unequality ||^(5)(icj) — ^(5i)(ia;)||2 < o"i holds. Moreover, G{Si) has no poles in 
C>o. 

Proof. First we show that G{Ss,^){s) > Q[Si){s) for all s G C>o except for a finite 

number. Let T be as in Lemma 3.11. We rewrite G{Ss,j){s) for 7 > ui and s £ C>o as 

giSs„){s) = GiTE-'^ ■ 5s„ • T-'){s) = (Ci 0) (TE-^isSs.-, - A..)T-i)"' r^i") + D 

-1 



c^i o)((^""-"" ^ +(.?-7^)(./ + ri?-TATT-i) rj+A 



= fr* n^ ( ( sEii — All 0\ ^ ^_2 „2\^„r , T^r-T /|T^^-l^ \ /^-B; 
where (*) holds with 



£3,1 = £s,^i + £3,1 - £3,^1 = £3,^1 + (CTI - 7 )^ ; 

A.^ = A,al-(^?-7^)^^• 
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We set I J.} J..? ) := TE "^ A^T ^ with appropriate dimensioned Mj. In addition, let 



._i _ (sEu - Au + {al - j^){sl + Ml) (af - ^^)M2 



V Ms Mi 

sia{M4)VJa{Eii,Aii). Then 

L^ := s^ii - 111 + {(yl - l^)isl + Ml) - {al - ^'^)M2{sI - M^y^M; 
is the Schur complement of 

TE-^is£s,, - As,,)T-' = ( -'' '^'\^, 

in (af — 7^)(s/ + M4). By Theorem 2.5 the matrix L-y is regular because {s£s,-y "-^3.7) is 
regular. The functions 7 1— )• L^ and thus 7 1— )• G{Ss,^){s) = CiL^^Bi + D are continuous 
in ((Ti,oo). Note that Lo-^ = sE'n — j4ii is regular. We obtain 

g{Ss,,){s) = CiL-^Bi + D ^^^ CiisEn - MiY^B^ + D = g{Si){s). 

Therefore for aU iw ^ aiM^) U a{Eii,An) 

\\g{s){iu)-g{Si){iu)\\^= lim \\g{s){ioj)-g{Ss,^,){iu)\\^< iim 7 = ^1. 

Now we show that ^(5*1) has no poles in C>o. By [Fra87, Section 2.3] we have ||G(s)||2 < 
||G||qq for all G £ RH^o and s G C>o. Thus we get the upper estimate 

\\g{Ss,.){s)h < mss^.){s) - g{s){s)\\, + \\g{s){s)\\, < 7 + \\g{s)\\^ . 

for all s £ C>o. Hence for all s £ C>o\(cr(-Eii,^ii) U 0"(M4)) we have 

\\g{Si){s)\\,= lim \\g{Ss,,){s)\\2< ^^Ji + \\Q{s)\U = cT, + \\g{s)\\^ . 

Thus ^(^i) has no poles in C>o. □ 

Theorem 3.13. The system Si given in Lemma 3.12 solves (APoo). 

Proof. We use the notations of the proof of Lemma 3.11. By Theorem 3.7 and Lemma 3.12 
it remains to prove that a{Eii, An) C C<oU{oo}. Let sq G C>o. We apply Lemma 2.15 

because sq is not a pole of g{Si). Let x = { ] satisfy = Cix and 



soT + Aj^r + Cj Ci^i -aiCjC2\ (xi 
-G1B2BJ -aiCjC2 \X2 



= (soE,, - All). = r°^ - ^- e" 7 ^' ' Trh),' \\1]- (3-2) 
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Hence 



^ = {Ci J=Cs,..,T2-M 1, (3.3) 






Due to 5;, ~ 5" G Sn,p,m "^^ have sq ^ (t(— A([) and the equation 

T 0\ /xi^ 



-, — ]^ I dj \ — , — — , — 1 I Jj 



follows. Thus we get xi = by the regularity of F. Now we deduce 

(3.2) ^ (TiC^d2X2 = ^ X2 = 

from the regularity of CjC2- Thus x = 0. Analogously every x G W^~ which fulfils 
= {soEii — Aii)~^ X and = BJx equals zero. Thus sq ^ a{Eu,Aii) and a{Eii, An) C 
C<oU{oo}. D 

The optimal solution given in the theorem above is still based on the balanced realization. 
Our second main contribution is the following result, which only needs one singular value 
decomposition. 

Theorem 3.14. Let S = {E,A,B,C,D) G S^^p^m- Suppose [£s,cr^,As,crj^) given in The- 
orem 3.1 is singular. Then there exist orthogonal U,V £ M"^" with 

UAs,^,V = Mji ^J^A ^^^ ^^g^i^^ Xu ■ (3.4) 

Every regular U,V G M"^" satisfying (3.4), fulfil 

The system S2 ■= (Eii,Aii,I3i,Ci,D] solves (APoo). 
Proof. Let T G M"^" be as in Lemma 3.11. Then we obtain 

{te-^£s^.,t-\te-^as^^,t-^) = ((^j^ °) , (^j^ ;;)) 
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with All G M("-'^)^("-'^) anda{Eu,Aii) C C<oU{oo}. In particular, An is regular. The 
matrices As,^^ and TE^~^As,aiT~^ have the same rank. By singular value decomposition 
there exist regular orthogonal U,V €z M*^^" such that 

UAs,.,V = r^Ji ^^A with regular An G m("-^)x("-'^) . 
Let U,V £ M"X" be regular and satisfy (3.4). We define 

For all s G C>o the matrices (si?ii — An) and (sE'n — j4ii) are regular. Thus we have 
fsEn-An ^Eu-AuV^^^^/^j^_,(sEn-An 0\ tv] [R^] 

y SE21 SE22 J V V 0/ / 

^^ (ue'^T-^) Tr"-'^ X {0}^] = R"-'^ X {0}^ . (3.5) 

We obtain the last equality because (*) holds for all s G C>o and in particular for s = 0. 
That implies E21 = 0, -£22 = and finally (T{En,An) = o"(-Eii, j4ii) which means 
5*2 G ^n-dpm- Furthermore, we have B2 = because 

Finally, we show that Q{S2) = Q{Si). Let s G p{En,An) and Yi := Ci(s£'ii — An)^^- 
We conclude with Y := (Yi O) G R^^" and X := (Xi X2) := YTE-~^U-^ that 

(ci a2)=(ci o)tv = yM^^-^^' ^\tv = x('^''-^^' sE,2-A,2 

In particular Xi = Ci{sEn — ^ii)^^- That implies 

g{Si){s) = CiisEn - An)-'Bi = Y^Bi = (^ Y2) {^^\ = (Yi Y2) TE'^U-' (^A 

= (Xi X2) (^^A = XiBi = diisEn - An)-'Bi = g(52)(s) . 

By Theorem 3.13 we get \\g{S) - g{S2)\\^ = ai. D 

Under the additional condition E = I we can use the less complex Schur decomposition 
to attain (3.4) in Theorem 3.14. 
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Lemma 3.15. In the situation of Theorem 3.14, where E = I, we obtain regular ma- 
trices U and V := U~^ , such that (3.4) holds, by applying the Schur decomposition to 

A,.i ■ 

Proof. Let T G M"^" be as in Lemma 3.11. Then we obtain 

with 111 G M("-^)^("-^) and cr(^ii,lii) C C<o U {oo}. In particular, In is regular. 
The matrices As,^^ and TAs,aiT~^ have the same eigenvalues counting multiplicities. By 
the Schur decomposition there exists a regular (in particular an orthogonal) U G M"^" 
such that 

where An G m("^'^)x("^^) is regular and a{A22) = {0}. This implies rank(Aii) = 
rank(2ii) = Tank{UAs,aiU~'^) and thus A22 = 0. D 

In summary, we obtain the following result for (APqc)- 

Theorem 3.16. Let S G ^npm- Then there exist S^ G S^ ^^ and S- = {E^, A^, B^, C_,0) G 



'^n-,p,m such that 5 ~ 5+ © 5_. We set cji := ^ /max a{E_Ds_E-^s_), then 

^ inf^ \\g{S)-g{S)\\oo = ai. 

We define the matrices 

^s_,.i := Os.E.dis.El - ap, £s_,., ■= EITZs,^,, Bs_,^, := EIDs_B_, 

As_,CT-^ '.= —A_TZs,a^ — C_Cs_,ai, Cs_,ai :=C-€s_E_. 

If {£s_,a^,As_,a^) is regular, then S+® {£s_,a^,As_,rr^,Bs_,a^,Cs_,ai,^) solves (AP^o). 
If {£s_,cTi,As_,ai) is singular, then there exist orthogonal matrices U,V & M"^" with 

UAs_,^,V = r^Ji "^^A and regular In . (3.6) 



Every regular U,V €z M"^" satisfying (3.6), fulfil 

'En Ei2\ jj^ _ (B 

The system S+ © (^n, ^n, i^i, Ci, 0) solves (AP^o) 



^^.-,.^=(7 7, f/^._,., = r;), cs_^.,v=(d, a 
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4 Application 

Algorithm 

We rewrite our results in an algorithmic manner. Note that the solution of (AP2) is 
necessary to compute the solution of (APqo) and thus no extra computation is needed. 

Input: S e S°^p^„ 

Output: Sq £ S^ p,„ which solves (APg) for q G {2, 00}. 

1. Decompose S into S- G §,^_ „ ,„ and 5+ G §^ „ „, as in Theorem 3.1 with the help 
of Theorem 2.3. 

2. The system 6*2 := 5*+ solves (AP2) . 



3. Determine €:^€s_EI, £3 := £'IDs_ and ai = Jma.xa{D^€) . 

4. Compute i? := D^£ - cr^/ , E-~e1r, B:=DB_, 

5. If {E^A) is regular (test with Corollary 3.8), then (APqo) is solved by So 

S+®{E,A,B,C,0). 

6. If {E,A) is singular, then determine [/, VeM"''" (SVD) such that 



UAV= (^" ^12 j ^i^j^ regular An. 



Compute Eii,Bi,Ci of UEV = (^^^ ^^A , UB = (^M , dv == (Ci C2) , then 
(APoo) is solved by Soc ■= S+ ® {En, An, Bi,di,0) . 

For the numerical issues of step 1 we refer to [KD92] and of step 3 to [Pen98] . 

Both approximation techniques does not increase the order of the system, i. e. Ug < 
n. For optimal i?/f2-approximation we even have 712 < n. Another advantage is the 
less computational complexity compared to i^ffoo-approximation. The behaviour of the 
transfer functions G{S) and ^(^2) at infinity is the same, i.e. Q{S){oo) = Q{S2){co). 
This does not hold for optimal i?i?oo-approximation in general. However, the suboptimal 
approximation S+ ® Ss_,-y (choosing 7 > cji in Theorem 3.7) has this property. Since 
£s_,-t is regular, we have 

g(S)(oo) = e(5+)(oo) = e(5+)(oo) + g(5s_,,)(oo) = G{S+ ® 5s_,,)(oo) . 

Examples 

The following benchmarks are performed in Matlab®. We test the building model 
(Figure 1) and the clamped beam model (Figure 2), see [CD02] for details. Both models 
are described by stable standard systems with one input and one output, i.e. m = p = 1. 
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We apply (shifted) Arnoldi method to reduce the systems. The parameters (interpolation 
point s G C U {00} and reduced order k G N) are based on the examples in [ASGOl]. 
The computed reduced systems are unstable as listed in Table 1. Thus we apply the 
algorithm presented above. Since the transfer functions of the original systems (and 
hence of the reduced systems) vanish at infinity an optimal i?//oo-approximation would 
result in huge relative errors at high frequencies. That is why we compute suboptimal 
approximations (i.e. 7 > ai) to match the original transfer functions at infinity. 



model 


order 


model reduction 


reduced 


unstable 


max. real part 






method 


order 


poles 


of unstable poles 


building 


48 


Arnoldi, s = 00 


31 


2 


?«42.4 






RH2 


29 


- 


- 






RHoo, 7 = 1.001- (71 


31 


- 


- 


clamped 


348 


Arnoldi, s = 0.1 


13 


1 


«1.5 


beam 




RH2 


12 


- 


- 






RHod = imi-ai 


13 


- 


- 



Table 1: Summary of the results 

The error between the original systems and the reduced systems of Arnoldi method does 
not significantly chance after applying our stabilization algorithm as depicted in Figure 1 
and Figure 2. 





Figure 1: Frequency responses of the (reduced) building model (left) and the error sys- 
tems (right) 
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Figure 2: Frequency responses of the (reduced) clamped beam model (left) and the error 
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